home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dtrexc.f < prev    next >
Text File  |  1996-07-19  |  10KB  |  347 lines
  1.       SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
  2.      $                   INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     March 31, 1993
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          COMPQ
  11.       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DTREXC reorders the real Schur factorization of a real matrix
  21. *  A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
  22. *  moved to row ILST.
  23. *
  24. *  The real Schur form T is reordered by an orthogonal similarity
  25. *  transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
  26. *  is updated by postmultiplying it with Z.
  27. *
  28. *  T must be in Schur canonical form (as returned by DHSEQR), that is,
  29. *  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  30. *  2-by-2 diagonal block has its diagonal elements equal and its
  31. *  off-diagonal elements of opposite sign.
  32. *
  33. *  Arguments
  34. *  =========
  35. *
  36. *  COMPQ   (input) CHARACTER*1
  37. *          = 'V':  update the matrix Q of Schur vectors;
  38. *          = 'N':  do not update Q.
  39. *
  40. *  N       (input) INTEGER
  41. *          The order of the matrix T. N >= 0.
  42. *
  43. *  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
  44. *          On entry, the upper quasi-triangular matrix T, in Schur
  45. *          Schur canonical form.
  46. *          On exit, the reordered upper quasi-triangular matrix, again
  47. *          in Schur canonical form.
  48. *
  49. *  LDT     (input) INTEGER
  50. *          The leading dimension of the array T. LDT >= max(1,N).
  51. *
  52. *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
  53. *          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
  54. *          On exit, if COMPQ = 'V', Q has been postmultiplied by the
  55. *          orthogonal transformation matrix Z which reorders T.
  56. *          If COMPQ = 'N', Q is not referenced.
  57. *
  58. *  LDQ     (input) INTEGER
  59. *          The leading dimension of the array Q.  LDQ >= max(1,N).
  60. *
  61. *  IFST    (input/output) INTEGER
  62. *  ILST    (input/output) INTEGER
  63. *          Specify the reordering of the diagonal blocks of T.
  64. *          The block with row index IFST is moved to row ILST, by a
  65. *          sequence of transpositions between adjacent blocks.
  66. *          On exit, if IFST pointed on entry to the second row of a
  67. *          2-by-2 block, it is changed to point to the first row; ILST
  68. *          always points to the first row of the block in its final
  69. *          position (which may differ from its input value by +1 or -1).
  70. *          1 <= IFST <= N; 1 <= ILST <= N.
  71. *
  72. *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
  73. *
  74. *  INFO    (output) INTEGER
  75. *          = 0:  successful exit
  76. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  77. *          = 1:  two adjacent blocks were too close to swap (the problem
  78. *                is very ill-conditioned); T may have been partially
  79. *                reordered, and ILST points to the first row of the
  80. *                current position of the block being moved.
  81. *
  82. *  =====================================================================
  83. *
  84. *     .. Parameters ..
  85.       DOUBLE PRECISION   ZERO
  86.       PARAMETER          ( ZERO = 0.0D+0 )
  87. *     ..
  88. *     .. Local Scalars ..
  89.       LOGICAL            WANTQ
  90.       INTEGER            HERE, NBF, NBL, NBNEXT
  91. *     ..
  92. *     .. External Functions ..
  93.       LOGICAL            LSAME
  94.       EXTERNAL           LSAME
  95. *     ..
  96. *     .. External Subroutines ..
  97.       EXTERNAL           DLAEXC, XERBLA
  98. *     ..
  99. *     .. Intrinsic Functions ..
  100.       INTRINSIC          MAX
  101. *     ..
  102. *     .. Executable Statements ..
  103. *
  104. *     Decode and test the input arguments.
  105. *
  106.       INFO = 0
  107.       WANTQ = LSAME( COMPQ, 'V' )
  108.       IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
  109.          INFO = -1
  110.       ELSE IF( N.LT.0 ) THEN
  111.          INFO = -2
  112.       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  113.          INFO = -4
  114.       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
  115.          INFO = -6
  116.       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
  117.          INFO = -7
  118.       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
  119.          INFO = -8
  120.       END IF
  121.       IF( INFO.NE.0 ) THEN
  122.          CALL XERBLA( 'DTREXC', -INFO )
  123.          RETURN
  124.       END IF
  125. *
  126. *     Quick return if possible
  127. *
  128.       IF( N.LE.1 )
  129.      $   RETURN
  130. *
  131. *     Determine the first row of specified block
  132. *     and find out it is 1 by 1 or 2 by 2.
  133. *
  134.       IF( IFST.GT.1 ) THEN
  135.          IF( T( IFST, IFST-1 ).NE.ZERO )
  136.      $      IFST = IFST - 1
  137.       END IF
  138.       NBF = 1
  139.       IF( IFST.LT.N ) THEN
  140.          IF( T( IFST+1, IFST ).NE.ZERO )
  141.      $      NBF = 2
  142.       END IF
  143. *
  144. *     Determine the first row of the final block
  145. *     and find out it is 1 by 1 or 2 by 2.
  146. *
  147.       IF( ILST.GT.1 ) THEN
  148.          IF( T( ILST, ILST-1 ).NE.ZERO )
  149.      $      ILST = ILST - 1
  150.       END IF
  151.       NBL = 1
  152.       IF( ILST.LT.N ) THEN
  153.          IF( T( ILST+1, ILST ).NE.ZERO )
  154.      $      NBL = 2
  155.       END IF
  156. *
  157.       IF( IFST.EQ.ILST )
  158.      $   RETURN
  159. *
  160.       IF( IFST.LT.ILST ) THEN
  161. *
  162. *        Update ILST
  163. *
  164.          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
  165.      $      ILST = ILST - 1
  166.          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
  167.      $      ILST = ILST + 1
  168. *
  169.          HERE = IFST
  170. *
  171.    10    CONTINUE
  172. *
  173. *        Swap block with next one below
  174. *
  175.          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
  176. *
  177. *           Current block either 1 by 1 or 2 by 2
  178. *
  179.             NBNEXT = 1
  180.             IF( HERE+NBF+1.LE.N ) THEN
  181.                IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
  182.      $            NBNEXT = 2
  183.             END IF
  184.             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
  185.      $                   WORK, INFO )
  186.             IF( INFO.NE.0 ) THEN
  187.                ILST = HERE
  188.                RETURN
  189.             END IF
  190.             HERE = HERE + NBNEXT
  191. *
  192. *           Test if 2 by 2 block breaks into two 1 by 1 blocks
  193. *
  194.             IF( NBF.EQ.2 ) THEN
  195.                IF( T( HERE+1, HERE ).EQ.ZERO )
  196.      $            NBF = 3
  197.             END IF
  198. *
  199.          ELSE
  200. *
  201. *           Current block consists of two 1 by 1 blocks each of which
  202. *           must be swapped individually
  203. *
  204.             NBNEXT = 1
  205.             IF( HERE+3.LE.N ) THEN
  206.                IF( T( HERE+3, HERE+2 ).NE.ZERO )
  207.      $            NBNEXT = 2
  208.             END IF
  209.             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
  210.      $                   WORK, INFO )
  211.             IF( INFO.NE.0 ) THEN
  212.                ILST = HERE
  213.                RETURN
  214.             END IF
  215.             IF( NBNEXT.EQ.1 ) THEN
  216. *
  217. *              Swap two 1 by 1 blocks, no problems possible
  218. *
  219.                CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
  220.      $                      WORK, INFO )
  221.                HERE = HERE + 1
  222.             ELSE
  223. *
  224. *              Recompute NBNEXT in case 2 by 2 split
  225. *
  226.                IF( T( HERE+2, HERE+1 ).EQ.ZERO )
  227.      $            NBNEXT = 1
  228.                IF( NBNEXT.EQ.2 ) THEN
  229. *
  230. *                 2 by 2 Block did not split
  231. *
  232.                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
  233.      $                         NBNEXT, WORK, INFO )
  234.                   IF( INFO.NE.0 ) THEN
  235.                      ILST = HERE
  236.                      RETURN
  237.                   END IF
  238.                   HERE = HERE + 2
  239.                ELSE
  240. *
  241. *                 2 by 2 Block did split
  242. *
  243.                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
  244.      $                         WORK, INFO )
  245.                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
  246.      $                         WORK, INFO )
  247.                   HERE = HERE + 2
  248.                END IF
  249.             END IF
  250.          END IF
  251.          IF( HERE.LT.ILST )
  252.      $      GO TO 10
  253. *
  254.       ELSE
  255. *
  256.          HERE = IFST
  257.    20    CONTINUE
  258. *
  259. *        Swap block with next one above
  260. *
  261.          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
  262. *
  263. *           Current block either 1 by 1 or 2 by 2
  264. *
  265.             NBNEXT = 1
  266.             IF( HERE.GE.3 ) THEN
  267.                IF( T( HERE-1, HERE-2 ).NE.ZERO )
  268.      $            NBNEXT = 2
  269.             END IF
  270.             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
  271.      $                   NBF, WORK, INFO )
  272.             IF( INFO.NE.0 ) THEN
  273.                ILST = HERE
  274.                RETURN
  275.             END IF
  276.             HERE = HERE - NBNEXT
  277. *
  278. *           Test if 2 by 2 block breaks into two 1 by 1 blocks
  279. *
  280.             IF( NBF.EQ.2 ) THEN
  281.                IF( T( HERE+1, HERE ).EQ.ZERO )
  282.      $            NBF = 3
  283.             END IF
  284. *
  285.          ELSE
  286. *
  287. *           Current block consists of two 1 by 1 blocks each of which
  288. *           must be swapped individually
  289. *
  290.             NBNEXT = 1
  291.             IF( HERE.GE.3 ) THEN
  292.                IF( T( HERE-1, HERE-2 ).NE.ZERO )
  293.      $            NBNEXT = 2
  294.             END IF
  295.             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
  296.      $                   1, WORK, INFO )
  297.             IF( INFO.NE.0 ) THEN
  298.                ILST = HERE
  299.                RETURN
  300.             END IF
  301.             IF( NBNEXT.EQ.1 ) THEN
  302. *
  303. *              Swap two 1 by 1 blocks, no problems possible
  304. *
  305.                CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
  306.      $                      WORK, INFO )
  307.                HERE = HERE - 1
  308.             ELSE
  309. *
  310. *              Recompute NBNEXT in case 2 by 2 split
  311. *
  312.                IF( T( HERE, HERE-1 ).EQ.ZERO )
  313.      $            NBNEXT = 1
  314.                IF( NBNEXT.EQ.2 ) THEN
  315. *
  316. *                 2 by 2 Block did not split
  317. *
  318.                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
  319.      $                         WORK, INFO )
  320.                   IF( INFO.NE.0 ) THEN
  321.                      ILST = HERE
  322.                      RETURN
  323.                   END IF
  324.                   HERE = HERE - 2
  325.                ELSE
  326. *
  327. *                 2 by 2 Block did split
  328. *
  329.                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
  330.      $                         WORK, INFO )
  331.                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
  332.      $                         WORK, INFO )
  333.                   HERE = HERE - 2
  334.                END IF
  335.             END IF
  336.          END IF
  337.          IF( HERE.GT.ILST )
  338.      $      GO TO 20
  339.       END IF
  340.       ILST = HERE
  341. *
  342.       RETURN
  343. *
  344. *     End of DTREXC
  345. *
  346.       END
  347.